# MEDIPS batch processing:
# Performing MEDIPS.methylationProfiling for pre-processed data

library(MEDIPS)
library(BSgenome.Hsapiens.UCSC.hg19)

getwd()
# [1] "/usr/local/data/Data/NCH/MBD-Seq/MEDIP-Analysis"

# get list of files in directory
# here we're loading the stored results of the previous batch run
files <- list.files(pattern="glioma.*RData$")

# performing the analysis for FANTOM4 promoters of chr7 
# ROI.file <- ("/home/skurscheid/Data/Riken/FANTOM4/Fantom4Promoters.chr7.txt")
ROI.file <- ("/home/skurscheid/Data/NCH/HOX/ACGHmurat/CGH_chr7BACs.bed")

#####################################
# beginning of loop
#####################################

#loop through files, process and write results to individual RData files 
for (i in files) {
	sample <- gsub(".RData","",i)
	print(sample)
	
	# sleep for between 30 to 60 seconds to allow I/O cool-down between parallel processes
	Sys.sleep(sample(30:60,5))
	
	# loading the RData file
	load(i)
	
	# writing MEDIPS data as WIG file for displaying in genome browser
	# first rms values as output
	#MEDIPS.exportWIG(file = paste(sample,".rms.WIG",sep=""), data = MEDIPS.SET, raw = FALSE, descr = paste(sample,"RMS values, genome-wide", sep=""))
	#system("gzip ")
	#MEDIPS.exportWIG(file = paste(sample,".rpm.WIG",sep=""), data = MEDIPS.SET, raw = TRUE, descr = paste(sample,"RPM values, genome-wide", sep=""))

	#chr7.fantom4.promoters.profiling <- MEDIPS.methylProfiling(data1 = MEDIPS.SET, ROI_file = ROI.file, math = mean, select = 2)
	#chr7.500bpwindows.profiling <- MEDIPS.methylProfiling(data1 = MEDIPS.SET, chr="chr7", frame_size = 500, step= 100, math = mean, select = 2)
	
	
	chr7.BACs.sample.profile <- MEDIPS.methylProfiling(data1 = MEDIPS.SET, ROI_file = ROI.file, math = mean, select = 2)
	
	# extraction of ams (Absolute Methylation Score) values, as these are corrected for CpG densities (Coupling Factors)
	# and therefore allow comparison/profiling across different ROIs with differing CpG densities
	if (which(files == i) == 1){
		chr7.BACs.ams.profile <- as.matrix(cbind(chr7.BACs.sample.profile$ams_A))
		rownames(chr7.BACs.ams.profile) <- rownames(chr7.BACs.sample.profile)
		colnames(chr7.BACs.ams.profile)[which(files == i)] <- sample
	} else {
		chr7.BACs.ams.profile <- cbind(chr7.BACs.ams.profile, as.matrix(chr7.BACs.sample.profile$ams_A))
		colnames(chr7.BACs.ams.profile)[which(files == i)] <- sample
	}
	
	# sleep for between 30 to 60 seconds to allow I/O cool-down between parallel processes
	Sys.sleep(sample(30:60,5))
	
	write.csv(chr7.BACs.sample.profile, file=paste(sample,"_chr7_CGHBACsProfile.csv", sep=""))
	
	#write.csv(chr7.fantom4.promoters.profiling, file=paste(sample,"_chr7_fantom4Promoters.MEDIPS.csv", sep=""))
	#write.csv(chr7.500bpwindows.profiling, file=paste(sample,"_chr7_500bpWindows.MEDIPS.csv", sep=""))
	
	gc()
}
	
# write collated AMS results to CSV
write.csv(chr7.BACs.ams.profile, file="chr7.BACs.ams.profile.csv")
	
###########################################
# end of loop
###########################################

